

% track diving vessels from movies:
warning('off','all')


%% set preferences:
deinterleave = 0;       %set to 1 if movie file is interleaved and want to use channel1; 2 to use channel 2; 0 if not
dofast = 1;             %set to 1 if doing fast tracking with just the radon transform; set to 0 for ellipse fitting (slower)
consistentROI = 0;      %set to 1 to just copy the ROIs from the first scan onto all the other scans (i.e. *001.tif gets the same ROIs as *002.tif and *003.tif)
spatialSmooth = true;   %true if want to do spatial filtering to emphasize edges
kernsz = 15;            %in microns
logsig = 2.5;             %in microns
gaussig = 1;                %in microns
temporalsmooth = 1;     %number of frames of smoothing to do in TIME; if this is 1 then no smoothing
%objtype = '16x_v2';
dostabilize = true;             %adds a image stabilization step if necessary
doDenoise = true;
if 0      %allows you to skip the selection if the selection has already been done
%% select diving ROIs; can do this separately if you want but they'll be saved
batchSelectDivingROIs(deinterleave,1);
%% apply ROI selections for technical replicates:
if consistentROI
    fn = dir('*1_selectedROIs.mat');
    for n = 1:length(fn)
        load(fn(n).name);
        %identify other .tifs with the same name
        ind = strfind(fn(n).name,'1_selected');
        fn2 = dir([fn(n).name(1:(ind-1)),'*.tif']);
        for k = 1:length(fn2)
            if ~strcmp(fn2(k).name,fn(n).name)
                savename = strrep(fn2(k).name,'.tif','_selectedROIs.mat');
                movname = fn2(k).name;
                save(savename,'roi','movname')
            end
        end
    end
end
end

%% go through the individual movies and track each ROI
fn = dir('*selectedROIs.mat');
posthocAveragingGlobal = 1;        % if you want to do additional time averaging (e.g. if you accidentally had averaging set to 1 instead of 3)
for n = 1:length(fn)
    clear movname roi results
    results = struct('params',[],'ci',[],'minoraxisFit',[],'radonEstimate',[],'moviename',[]);
    load(fn(n).name);
    if ~exist('movname','var')
        movname = strrep(fn(n).name,'_selectedROIs.mat','.tif');
    end
    savename = strrep(movname,'.tif','_divingTracking6.mat');
    if isfile(savename)    %this statement prevents re-tracking already tracked files....
        continue;
    end
    meta = getSImetadata(movname);
    %pix2mic = 1000/meta.Height/meta.zoom;
    if meta.frameaveraging == 1 && round(meta.framerate,0) == 30
        posthocAveraging = 3;
    else
        posthocAveraging = posthocAveragingGlobal;
    end
    disp(['using settings for this objective: ',objtype]);
    pix2mic = pixel2micron_2p(meta.zoom,objtype);
    kernszPix = round(kernsz/pix2mic);
    logsigPix = logsig/pix2mic;
    gaussigPix = gaussig/pix2mic;
    [X,Y] = meshgrid(-kernszPix:kernszPix, -kernszPix:kernszPix);
    kern = (1/(pi*logsigPix^4))*...
        (1-.5*((X.^2 + Y.^2)/(logsigPix^2))).*...
        exp((-(X.^2 + Y.^2))/(2*logsigPix^2));
    kern = kern./sum(kern(:));

    if doDenoise  %this will run a denoising algorithm that seems to work reasonably well, but requires loading in the movie differently:
        fnRawImg = fullfile(cd,movname);
        channels = struct('blood_plasma',1);
        load('CHIPS_25xCalibration.mat');
        rawImg = SCIM_Tif(fnRawImg,channels,calibration);
        rawImg.denoise();
        mov = double(rawImg.rawdata);
        mov = reshape(mov,meta.Height,meta.Width,[]);
    end
    switch deinterleave
        case 1            
            if length(meta.savedchannels) < 2 & ~doDenoise
                mov = loadTiffStack2_SI(movname);
            elseif doDenoise
                nf = meta.totalframes/meta.frameaveraging*2;
                mov = mov(:,:,1:2:nf);
            else
                nf = meta.totalframes/meta.frameaveraging*2;
                mov = loadTiffStack2_SI(movname,1:2:nf);
            end
        case 0
            if ~doDenoise
                mov = loadTiffStack2_SI(movname);
            end
        case 2
            if length(meta.savedchannels) < 2 & ~doDenoise
                mov = loadTiffStack2_SI(movname);
            elseif doDenoise
                nf = meta.totalframes/meta.frameaveraging*2;
                mov = mov(:,:,2:2:nf);
            else
                nf = meta.totalframes/meta.frameaveraging*2;
                mov = loadTiffStack2_SI(movname,2:2:nf);
            end
    end
    if posthocAveraging > 1
        mov = movmean(mov,posthocAveraging,3);
        mov = mov(:,:,(ceil(posthocAveraging/2):posthocAveraging:end));
        fps = meta.framerate/posthocAveraging;
    end
    if dostabilize
        refframe = 1;
        mov = driftCorrectWidefield(mov,refframe);
    end
    roiPosition = zeros(size(roi,1),2);
    for k = 1:size(roi,1)
        curmov = mov(roi(k,1):roi(k,2),roi(k,3):roi(k,4),:);
        if temporalsmooth > 1
            curmov = movmedian(curmov,temporalsmooth,3);
        end
        if spatialSmooth
            for framecount = 1:size(curmov,3)
                temp = curmov(:,:,framecount);
                temp = medfilt2(temp);
                %temp = conv2(temp,kern,'same');
                temp = imgaussfilt(temp,gaussigPix);
                temp(temp<0) = 0;
                curmov(:,:,framecount) = temp;
            end
        end
        %[params,ci,minoraxis,radonEstimate] = fitEllipseFilledVessel(curmov,dofast);
        
        [params,ci,minoraxis,radonEstimate] = fitEllipseFilledVessel_v3_0(curmov,dofast,pix2mic);
        roiPosition(k,:) = round([roi(k,1)+roi(k,2), roi(k,3)+roi(k,4)]/2);         %[row,col] position of the center of the ROI in the movie
        results(k).params = params;
        results(k).ci = ci;
        results(k).minoraxisFit = minoraxis;
        results(k).radonEstimate = radonEstimate;
        results(k).moviename = movname;
        results(k).objtype = objtype;
    end
    
    
    %% extract position information:
    %meta = getSImetadata(movname);
    rootposition = repmat(meta.position,size(roiPosition,1),1);
    %pix2mic = 1000/512/meta.zoom;
    
    %switch ROI position to x,y from row,col
    roiPosition = [roiPosition(:,2),roiPosition(:,1)];
    
    offset = [meta.Height/2,meta.Width/2]*pix2mic;
    roiPosition = roiPosition*pix2mic-offset;     %convert pixel position to micron position
    roiPosition = [roiPosition + rootposition(:,1:2),rootposition(:,3)];
    for k = 1:length(results)
        results(k).roiPosition = roiPosition(k,:);
    end
    save(savename,'results')
    disp(['finished ',num2str(n/length(fn))])
end
    
warning('on','all')
